Consider the non-homogeneous vector wave equation,
\begin{align}
	\nabla^2\vec{A}(\vec{r})-\gamma^2\vec{A}(\vec{r})=\nabla^2\vec{A}(\vec{r})+k^2\vec{A}(\vec{r})=\vec{J}(\vec{r})\label{eqn:nonhomowaveeqn}
\end{align}
where $\vec{A}$ is a vector function of spatial coordinates $\vec{r}$.  The vector function $\vec{J}$ is a source or driving function and $\gamma^2=-k^2$ is a complex propagation constant determined by the medium in which the waves are propagating.  This wave equation is identical to (\ref{eqn:waveeqnsoln}) with $\vec{J}=0$, which results in a homogeneous wave equation.  In order to keep the discussion general, the source $\vec{J}$ was added, which will become more relevant when discussing vector potential theory.

The vector functions $\vec{A}$ and $\vec{J}$ are separated into
transverse $(\bot)$ and normal $(n)$ components as follows.
\begin{align}
    \vec{A}&=\vec{A}_\bot+A_n\hat{n}\label{eqn:sepvec}\\
    \vec{J}&=\vec{J}_\bot+A_n\hat{n}\label{eqn:sepvecJ}
\end{align}
The operator $\vec{\nabla}$ and the Laplacian operator are also separated as follows,
\begin{align}
    \vec{\nabla}&=\vec{\nabla}_\bot+\frac{\partial}{\partial n}\hat{n}\label{eqn:sepnabla}\\
    \nabla^2&=\vec{\nabla}\cdot\vec{\nabla}=\nabla_\bot^2+\frac{\partial^2}{{\partial{n}}^2}\label{eqn:sepnablasqr}
\end{align}
Expanding the wave equation (\ref{eqn:nonhomowaveeqn}) into transverse
and normal components yields,
\begin{equation}\label{eqn:vecwavetnexp1}
\begin{split}
\left(\nabla_\bot^2+\frac{\partial^2}{{\partial{n}}^2}\right)(\vec{A}_\bot+A_n\hat{n})+(\vec{J}_\bot+J_n\hat{n})&=\gamma^2(\vec{A}_\bot+A_n\hat{n})\\
&=-k^2(\vec{A}_\bot+A_n\hat{n})
\end{split}
\end{equation}
Expanding and equating components gives the following two wave
equations,
\begin{equation}\label{eqn:vecwaveperp}
    \left(\nabla^2_\bot+\frac{\partial^2}{{\partial n}^2}\right)\vec{A_\bot}+\vec{J}_\bot=\gamma^2\vec{A}_\bot=-k^2\vec{A}_\bot
\end{equation}
\begin{equation}\label{eqn:vecwavenorm}
    \left(\nabla^2_\bot+\frac{\partial^2}{{\partial n}^2}\right)A_n+J_n=\gamma^2A_n=-k^2A_n
\end{equation}
Equation (\ref{eqn:vecwaveperp}) is a two dimensional vector wave
equation and (\ref{eqn:vecwavenorm}) is a scalar wave equation.
The \emph{separation of variables} technique will now be utilized
to solve equation (\ref{eqn:vecwavenorm}).  Let $A_n(\vec{r})$ and $J_n(\vec{r})$,
which is a function of three independent variables, be separated
into the product of two functions, with one function being
dependent on the transverse variables and the other function being
dependent on the normal variable.  Let the transverse function
with its independent variable be denoted as $a_n(\vec{r_\bot})$ and $j_n(\vec{r_\bot})$ and
the normal function with its independent variable as $N(n)$. In
equation form this is expressed as,
\begin{align}
    A_n(\vec{r})&=a_n(\vec{r}_\bot)N(n)\label{eqn:sepfcns}\\
    J_n(\vec{r})&=j_n(\vec{r}_\bot)N(n)\label{eqn:sepfcnsjn}
\end{align}
Substituting (\ref{eqn:sepfcns}) and (\ref{eqn:sepfcnsjn}) into (\ref{eqn:vecwavenorm})
yields,
\begin{equation}\label{eqn:vecwavenorm2}
\begin{split}
N(n)\nabla^2_\bot a_n(\vec{r_\bot})+a_n(\vec{r_\bot})\frac{d^2N(n)}{{d n}^2}+j_n(\vec{r}_\bot)N(n)&=\gamma^2a_n(\vec{r_\bot})N(n)\\
&=-k^2a_n(\vec{r_\bot})N(n)
\end{split}
\end{equation}
Note that the transverse operator $\nabla_\bot^2$ only operates on
the transverse function $a_n(\vec{r_\bot})$. Likewise, the normal
operator $\frac{d^2}{{d n}^2}$ only operates on the normal
function $N(n)$.  The partial derivative was changed to an
ordinary derivative, since it is operating only on a single
variable. Dividing both sides of (\ref{eqn:vecwavenorm2}) by
$a_n(\vec{r_\bot})N(n)$ yields,
\begin{equation}\label{eqn:vecwavenorm3}
\left(\frac{\nabla^2_\bot a_n(\vec{r_\bot})}{a_n(\vec{r_\bot})}+\frac{j_n}{a_n}\right)+\left(\frac{1}{N(n)}\frac{d^2N(n)}{{d n}^2}\right)=\gamma^2=-k^2
\end{equation}
The left hand side of (\ref{eqn:vecwavenorm3}) is a summation of
two terms, each of which are dependent on distinct variables,
which sum to equal a constant.  The only way for this to be true,
is that each term equal a constant. Therefore,
\begin{equation}\label{eqn:scalarwaveperp0}
    \frac{\nabla^2_\bot a_n(\vec{r_\bot})}{a_n(\vec{r_\bot})}=\gamma_\bot^2=-k_\bot^2
\end{equation}
\begin{equation}\label{eqn:scalarwavenorm0}
    \frac{1}{N(n)}\frac{d^2N(n)}{{d n}^2}=\gamma_n^2=-k_n^2
\end{equation}
This leads to the following differential equations.
\begin{equation}\label{eqn:scalarwaveperp}
    \nabla^2_\bot a_n(\vec{r_\bot})+j_n=\gamma_\bot^2a_n(\vec{r_\bot})=-k_\bot^2a_n(\vec{r_\bot})
\end{equation}
\begin{equation}\label{eqn:scalarwavenorm}
    \frac{d^2N(n)}{{d n}^2}=\gamma_n^2N(n)=-k_n^2N(n)
\end{equation}
Because each of the terms equal a constant and sum to equal a
constant the following equation is true, known as the
\emph{constraint equation}.
\begin{equation}\label{eqn:constrainteqn2d}
    \boxed{\gamma^2=\gamma_n^2+\gamma_\bot^2=-k^2=-k_n^2-k_\bot^2}
\end{equation}

Now consider the following important assumption that the transverse vector fields $\vec{A}_\bot(\vec{r})$ and $\vec{J}_\bot(\vec{r})$ can be separated into the product of a transverse terms $\vec{a}(\vec{r}_\bot)$ and $\vec{j}(\vec{r}_\bot)$ and a normal term $N(n)$, or
\begin{align}
	\vec{A}_\bot(\vec{r})=\vec{a}(\vec{r}_\bot)N(n)\label{eqn:sepfcns2}\\
	\vec{J}_\bot(\vec{r})=\vec{j}(\vec{r}_\bot)N(n)\label{eqn:sepfcns2jn}
\end{align}
in which $N(n)$ was defined previously in (\ref{eqn:sepfcns}). Substituting (\ref{eqn:sepfcns2}) and (\ref{eqn:constrainteqn2d}) into (\ref{eqn:vecwaveperp}) gives,
\begin{align}
\begin{split} \left(\nabla^2_\bot+\frac{\partial^2}{{\partial{n}}^2}\right)\vec{a}(\vec{r}_\bot)N(n)+\vec{j}(\vec{r}_\bot)N(n)&=(\gamma_n^2+\gamma_\bot^2)\vec{a}(\vec{r}_\bot)N(n)\\
&=(-k_n^2-k_\bot^2)\vec{a}(\vec{r}_\bot)N(n)
\end{split}
\end{align}
and then expanding and factoring gives,
\begin{align} 
N(n)\left(\nabla^2_\bot\vec{a}(\vec{r}_\bot)-\gamma_\bot^2\vec{a}(\vec{r}_\bot)\right)+\vec{a}(\vec{r}_\bot)\left(\frac{\partial^2N(n)}{{\partial{n}}^2}-\gamma^2N(n)\right)+\vec{j}(\vec{r}_\bot)N(n)=0\label{eqn:vectorwaveperp1}
\end{align}
Finally, substituting (\ref{eqn:scalarwavenorm}) into (\ref{eqn:vectorwaveperp1}) and simplifying yields a transverse wave equation that is only a function of the transverse independent variables,
\begin{align} 
\nabla^2_\bot\vec{a}(\vec{r}_\bot)-\gamma_\bot^2\vec{a}(\vec{r}_\bot)=\nabla^2_\bot\vec{a}(\vec{r}_\bot)+k_\bot^2\vec{a}(\vec{r}_\bot)=\vec{j}(\vec{r}_\bot)\label{eqn:vectorwaveperp}
\end{align}

In order to arrive at these results, the field variation in the normal direction was factored out of the vector fields (both the normal component and the transverse components), and then the vector field was broken down into normal and transverse components.  Therefore, the fields and sources can be expressed as,
\begin{align}
	\vec{A}(\vec{r}) = \left[\vec{a}(\vec{r}_\bot)+a_n(\vec{r}_\bot)\hat{n}\right]N(n)\label{eqn:transmissionlinefield}\\
	\vec{J}(\vec{r}) = \left[\vec{j}(\vec{r}_\bot)+j_n(\vec{r}_\bot)\hat{n}\right]N(n)
\end{align}
The constant $\gamma_\bot^2=-k_\bot^2$ can be found by applying boundary conditions and solving (\ref{eqn:scalarwaveperp}) or (\ref{eqn:vectorwaveperp}) and then the constant $\gamma_n^2=-k_n^2$ can then be found from the constraint equation (\ref{eqn:constrainteqn2d}).